!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! 
!! Turbulence closures
!!
!! \n
!!
!! This module provides a collection of turbulence closures for the
!! (weakly) compressible Navier-Stokes equations. The turbulence
!! fluxes are computed pointwise, independently of the grid and the
!! finite element basis. Gradients of the relevant fields must be
!! computed elsewhere.
!!
!! The turbulent fluxes are collected in a three dimensional array
!! \f${\bf F}_{turb}\f$ as follows:
!! \f{displaymath}{
!!  \begin{array}{ccl}
!!  {\bf F}_{turb}(\,:\,,\phantom{2:}1\phantom{+d},l) & = &
!!  \underline{\mathcal{F}}^d_E(x_l)
!!  \\[1mm]
!!  {\bf F}_{turb}(\,:\,,2:1+d,l) & = &
!!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}}(x_l)
!!  \end{array}
!! \f}
!! where \f$x_l\f$, \f$l=1,\ldots,N_{points}\f$ are the coordinates of
!! the points where the flux is required.
!!
!! \note The functions provide also some optional arguments which can
!! be used to recover intermediate values used in the computation of
!! the turbulent fluxes.
!!
!! \section uniform_expansion Uniform expansion
!!
!! It is useful to derive the relation between the symmetric velocity
!! gradient and the velocity divergence for the case of a uniform
!! expansion, since sometimes the dissipative fluxes are defined so
!! that they are zero for such a flow. Let us denote by \f$d\f$ the
!! spatial dimension and let us consider a flow characterized by:
!! \f{displaymath}{
!!  \nabla\underline{u} = \eta\,\mathcal{I}
!! \f}
!! for some expansion velocity \f$\eta\f$. We have that
!! \f{displaymath}{
!!  \nabla\underline{u}+\nabla\underline{u}^T = 2\eta\, \mathcal{I} =
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
!! \f}
!! Hence, to guarantee that a uniform expansion does not generate any
!! dissipative stress, we can express the dissipative fluxes as a
!! function of the deviatoric deformation velocity
!! \f{displaymath}{
!!  \mathcal{S}_D = \nabla\underline{u}+\nabla\underline{u}^T -
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
!! \f}
!! 
!! \subsection Example
!!
!! Consider the flow
!! \f{displaymath}{
!!  u_1 = \eta x_1,\qquad u_2 = \eta x_2, \qquad u_3 = 0.
!! \f}
!! If this flow is seen as a two-dimensional flow, \f$\mathcal{S}_D\f$
!! vanishes, however if it is regarded as a three-dimensional flow we
!! have
!! \f{displaymath}{
!!  \mathcal{S}_D = \frac{\eta}{3}\left[ 
!!   \begin{array}{ccc}
!!     2 & 0 & 0 \\
!!     0 & 2 & 0 \\
!!     0 & 0 & -4 \\
!!   \end{array}
!!   \right].
!! \f}
!!
!<----------------------------------------------------------------------
module mod_turb_flux

!-----------------------------------------------------------------------

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   t_phc, phc, coeff_visc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_turb_flux_constructor, &
   mod_turb_flux_destructor,  &
   mod_turb_flux_initialized, &
   i_turb_mod, t_turb_input, t_turb_diags, &
   viscous_flux, smagrich_flux

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> turbulence model input
 type t_turb_input
   !> \f$\nabla\underline{u}+\nabla\underline{u}^T\f$
   real(wp), allocatable :: gsuu(:,:,:)
   !> temperature gradient
   real(wp), allocatable :: gtt(:,:)
   !> potential temperature gradient
   real(wp), allocatable :: gth(:,:)
   !> tracer gradients (second index for tracer \c it)
   real(wp), allocatable :: gcc(:,:,:)
   !> \f$\displaystyle \frac{1}{|K|}\int_K \sum_{ij}\left(
   !! \nabla\underline{u}+\nabla\underline{u}^T
   !! \right)_{ij}^2\,d\underline{x}\f$
   real(wp) :: s2_e
   !> element average vertical gradient of potential temperature
   real(wp) :: gzth_e
 end type t_turb_input

 !> turbulence diagnostics
 !!
 !! Type which groups all the possible diagnostics returned by the
 !! various turbulence models
 !<
 type t_turb_diags
   real(wp), allocatable :: nu(:)    !< turbulent viscosity
   real(wp), allocatable :: ri(:)    !< Richardson number
   real(wp), allocatable :: s_abs(:) !< \f$ |\mathcal{S}| \f$
   !> turbulent dissipation \f$
   !! - \underline{\underline{\mathcal{F}}}^d_{\underline{U}} :
   !! \mathcal{S} \f$
   !<
   real(wp), allocatable :: diss(:)
   !> mean deformation rate
   real(wp), allocatable :: s2_e(:)
 end type t_turb_diags

 !> This is the general interface for all the turbulence models
 abstract interface
  pure subroutine i_turb_mod(fem , d,h,x,rho,p,uu , ti , turb_diags)
   import :: wp, t_turb_input, t_turb_diags
   !> spatial dimension
   integer, intent(in) :: d
   !> element size
   real(wp), intent(in) :: h
   !> positions at which the fluxes are required
   real(wp), intent(in) :: x(:,:)
   !> density
   real(wp), intent(in) :: rho(:)
   !> pressure
   real(wp), intent(in) :: p(:)
   !> velocity
   real(wp), intent(in) :: uu(:,:)
   !> additional input
   type(t_turb_input), intent(in) :: ti
   !> energy and momentum dissipative fluxes
   real(wp), intent(out) :: fem(:,:,:)
   !> Turbulence diagnostics
   type(t_turb_diags), intent(out), optional :: turb_diags
  end subroutine i_turb_mod
 end interface

! Module variables

 ! public members
 logical, protected ::               &
   mod_turb_flux_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_turb_flux'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_turb_flux_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
(mod_dgcomp_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_turb_flux_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_turb_flux_initialized = .true.
 end subroutine mod_turb_flux_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_turb_flux_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_turb_flux_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_turb_flux_initialized = .false.
 end subroutine mod_turb_flux_destructor

!-----------------------------------------------------------------------

 !> Viscous flux, prescribed scalar viscosity
 !!
 !! The energy and momentum diffusive fluxes are
 !! \f{displaymath}{
 !!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
 !!  -\rho\nu\left[ \nabla\underline{u} + \nabla\underline{u}^T +
 !!  \lambda\nabla\cdot\underline{u}\,\mathcal{I}\right], \qquad
 !!  \underline{\mathcal{F}}^d_E = -\rho\frac{\nu c_p}{Pr}\nabla T +
 !!  \underline{u}\cdot
 !!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}}.
 !! \f}
 !<
 pure subroutine viscous_flux(fem , d,h,x,rho,p,uu , ti , turb_diags)
  integer, intent(in) :: d
  real(wp), intent(in) :: h !< unused
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in) :: rho(:), p(:) !< p unused
  real(wp), intent(in) :: uu(:,:)
  type(t_turb_input), intent(in) :: ti
  real(wp), intent(out) :: fem(:,:,:)
  type(t_turb_diags), intent(out), optional :: turb_diags

  integer :: l, id
  real(wp) :: nu(2,size(x,2)), ldivu, sym(d,d)
 
   ! Compiler bug: this should not be necessary, however see
   ! http://software.intel.com/en-us/forums/showthread.php?t=82775&o=a&s=lr
   ! and bug DPD200169598 for ifort.
   if(present(turb_diags)) then
     if(allocated(turb_diags%s_abs)) deallocate(turb_diags%s_abs)
     if(allocated(turb_diags%diss )) deallocate(turb_diags%diss )
   endif
   ! end of compiler bug

   if(present(turb_diags)) then
     allocate(turb_diags%s_abs(size(x,2)))
     allocate(turb_diags%diss (size(x,2)))
   endif

   ! problem coefficients
   nu = coeff_visc( x , uu )

   do l=1,size(x,2)
     ! momentum flux
     ldivu = (-2.0_wp/3.0_wp)*0.5_wp*tr(ti%gsuu(:,:,l)) ! lambda*div(u)
     sym = ti%gsuu(:,:,l) ! symmetric gradient
     do id=1,d
       sym(id,id) = sym(id,id) + ldivu
     enddo
     fem(:,2:1+d,l) = -rho(l)*nu(2,l)*sym

     ! energy flux
     fem(:,  1  ,l) = -rho(l)*nu(1,l)*phc%cp*ti%gtt(:,l) &
                      + matmul( uu(:,l) , fem(:,2:1+d,l) )

     ! tracers
     fem(:,2+d: ,l) = -rho(l)*nu(1,l)*ti%gcc(:,:,l)

     if(present(turb_diags)) then
       turb_diags%s_abs(l) = sqrt(sum(sym**2))
       turb_diags%diss (l) = -sum(fem(:,2:1+d,l)*sym) ! use sym^T=sym
     endif
   enddo

 end subroutine viscous_flux

!-----------------------------------------------------------------------

 !> Smagorinsky-Richardson isotropic turbulent viscosity
 !!
 !! In this subroutine we consider the model discussed in <a
 !! href="http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1962.tb00128.x/pdf">
 !! [Lilly, 1962]</a> and <a
 !! href="http://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%281983%29111%3C2341%3AACMFTS%3E2.0.CO%3B2">
 !! [Durran, Klemp, 1983]</a>. In practice, this amounts to a
 !! Smagorinsky eddy viscosity corrected according to the local
 !! Richardson number.
 !!
 !! We set
 !! \f{displaymath}{
 !!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
 !!  -\rho\nu_U\mathcal{S}, \qquad
 !!  \underline{\mathcal{F}}^d_E = -c_p\rho\nu_E \nabla T.
 !! \f}
 !! with
 !! \f{displaymath}{
 !!  \mathcal{S} = \nabla\underline{u} + \nabla\underline{u}^T -
 !!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
 !! \f}
 !! The turbulent viscosities are
 !! \f{displaymath}{
 !!  \nu_U = k^2 \Delta x\,\Delta z\, \frac{|\mathcal{S}|}{\sqrt{2}}
 !!  \sqrt{\max\left(1-\frac{Ri}{Pr}\,,0\right)}, \qquad \nu_E =
 !!  Pr^{-1}\nu_U,
 !! \f}
 !! where \f$k=0.21\f$, \f$Pr=\frac{1}{3}\f$ and
 !! \f{displaymath}{
 !!  |\mathcal{S}| = \sqrt{\sum_{i,j}\mathcal{S}^2_{ij}} , \qquad Ri =
 !!  2 \frac{\frac{g}{\theta} \frac{\partial\theta}{\partial
 !!  z}}{|S|^2}.
 !! \f}
 !<
 pure subroutine smagrich_flux(fem , d,h,x,rho,p,uu , ti , turb_diags)
  integer, intent(in) :: d
  real(wp), intent(in) :: h
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in) :: rho(:), p(:)
  real(wp), intent(in) :: uu(:,:)
  type(t_turb_input), intent(in) :: ti
  real(wp), intent(out) :: fem(:,:,:)
  type(t_turb_diags), intent(out), optional :: turb_diags

  logical, parameter :: local_ri = .false.

  real(wp), parameter :: &
    k  = 0.21_wp,       &
    pr = 1.0_wp/3.0_wp
  integer :: l, id
  real(wp) :: theta(size(x,2)), ldivu, sym(d,d), ri, mu

   ! Compiler bug: this should not be necessary, however see
   ! http://software.intel.com/en-us/forums/showthread.php?t=82775&o=a&s=lr
   ! and bug DPD200169598 for ifort.
   if(present(turb_diags)) then
     if(allocated(turb_diags%s_abs)) deallocate(turb_diags%s_abs)
     if(allocated(turb_diags%ri   )) deallocate(turb_diags%ri   )
     if(allocated(turb_diags%nu   )) deallocate(turb_diags%nu   )
     if(allocated(turb_diags%diss )) deallocate(turb_diags%diss )
     if(allocated(turb_diags%s2_e )) deallocate(turb_diags%s2_e )
   endif
   ! end of compiler bug

   if(present(turb_diags)) then
     allocate(turb_diags%s_abs(size(x,2)))
     allocate(turb_diags%ri   (size(x,2)))
     allocate(turb_diags%nu   (size(x,2)))
     allocate(turb_diags%diss (size(x,2)))
     allocate(turb_diags%s2_e (size(x,2)))
   endif
 
   theta = p/(rho*phc%rgas) * (phc%p_s/p)**phc%kappa
   do l=1,size(x,2)
     ! momentum flux
     ldivu = (-2.0_wp/real(d,wp)) * 0.5_wp*tr(ti%gsuu(:,:,l))
     sym = ti%gsuu(:,:,l) ! symmetric gradient
     do id=1,d
        sym(id,id) = sym(id,id) + ldivu
     enddo

     if(local_ri) then
       ! Richardson number multiplied by |sym|^2 to avoid 0 division
       ri = 2.0_wp * phc%gravity*ti%gth(d,l)/theta(l)
       mu = rho(l) * (k**2/sqrt(2.0_wp)) * h**2  &
             * sqrt(max( sum(sym**2)-ri/pr , 0.0_wp ))
     else
       ! all the terms including derivatives are averaged over the
       ! gradient to compensate the numerical instability of 
       ri = 2.0_wp * phc%gravity*ti%gzth_e/theta(l)
       mu = rho(l) * (k**2/sqrt(2.0_wp)) * h**2  &
             * sqrt(max( ti%s2_e-ri/pr , 0.0_wp ))
     endif

     fem(:,2:1+d,l) = -mu * sym

     ! energy flux
     fem(:,  1  ,l) = -phc%cp * mu/pr * ti%gtt(:,l)

     ! tracers
     fem(:,2+d: ,l) = -mu/pr * ti%gcc(:,:,l)

     if(present(turb_diags)) then
       turb_diags%s_abs(l) = sqrt(sum(sym**2))
       turb_diags%ri   (l) = ri ! more precisely:  |S|^2 * ri = 2*N^2
       turb_diags%nu   (l) = mu/rho(l)
       turb_diags%diss (l) = -sum(fem(:,2:1+d,l)*sym) ! use sym^T = sym
       turb_diags%s2_e (l) = sqrt(ti%s2_e)
     endif
   enddo

 end subroutine smagrich_flux

!-----------------------------------------------------------------------

 !> Anisotropic model.
!  !!
!  !! In this subroutine we follow the general framework of <a
!  !! href="http://www.sciencedirect.com/science?_ob=MImg&_imagekey=B6TYJ-4B82P8M-3M-1&_cdi=5620&_user=2620285&_pii=S0898122103900149&_origin=gateway&_coverDate=08%2F31%2F2003&_sk=999539995&view=c&wchp=dGLzVlz-zSkWb&md5=bd04c1f14bad52d9b4db33591284d7ee&ie=/sdarticle.pdf">
!  !! [Abb&agrave;, Cercignani, Valdettaro, 2003]</a> to prescribe a
!  !! turbulent diffusion in the vertical direction as discussed in <a
!  !! href="http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1962.tb00128.x/pdf">
!  !! [Lilly, 1962]</a>. In practice, this amounts to a Smagorinsky eddy
!  !! viscosity corrected according to the Richardson number and
!  !! restricted to the vertical direction.
!  !!
!  !! Let us first consider the momentum flux and let us write, for
!  !! simplicity, \f$\tau_{ij} = \left(
!  !! \underline{\underline{\mathcal{F}}}^d_{\underline{U}}
!  !! \right)_{ij}\f$. The idea is to write
!  !! \f{displaymath}{
!  !!  \tau_{ij} = -2 \tilde{\Delta}^2|S|B_{ijrs} S_{rs},
!  !! \f}
!  !! where \f$S = \nabla\underline{u} + \nabla\underline{u}^T\f$ and we
!  !! have kept the isotropic component. The diffusion tensor is now
!  !! assumed to have the following form:
!  !! \f{displaymath}{
!  !!  B_{ijrs} = \sum_{\alpha,\beta} C_{\alpha\beta}\, a_{i\alpha}\,
!  !!  a_{j\beta}\, a_{r\alpha}\, a_{s\beta}
!  !! \f}
!  !! where \f$a\f$ denotes the normal matrix whose columns are \f$d\f$
!  !! orthonormal vectors \f$\underline{a}_\alpha\f$ and
!  !! \f$C_{\alpha\beta}\f$ is a symmetric tensor. If we now consider a
!  !! reference system aligned with the unit vectors
!  !! \f$\underline{a}_\alpha\f$, so that \f$a_{ij}=\delta_{ij}\f$, we
!  !! obtain
!  !! \f{displaymath}{
!  !!  \tau_{ij} = -2 \tilde{\Delta}^2|S|C_{ij} S_{ij} \qquad {\rm
!  !!  (no\,\,sum)}.
!  !! \f}
!  !! To determine the values \f$C_{ij}\f$, let us consider the
!  !! hypothesis that the momentum diffusion only takes place along the
!  !! vertical direction, which we denote by \f$z\f$, i.e.
!  !! \f{displaymath}{
!  !!  \tau_{\mu\nu} = \delta_{z\}
!  !! \f}
 !<
! pure subroutine aniso_flux(fem , d,x,rho,uu,guu,gtt)
!  integer, intent(in) :: d
!  real(wp), intent(in) :: x(:,:)
!  real(wp), intent(in) :: rho(:)
!  real(wp), intent(in) :: uu(:,:), guu(:,:,:)
!  real(wp), intent(in) :: gtt(:,:)
!  real(wp), intent(out) :: fem(:,:,:)
!
!  integer :: l, id
!  real(wp) :: nu(2,size(x,2)), ldivu, sym(d,d)
! 
!   ! problem coefficients
!   nu = coeff_visc( x , uu )
!
!   do l=1,size(x,2)
!     ! momentum flux
!     ldivu = (-2.0_wp/3.0_wp) * tr(guu(:,:,l)) ! lambda*div(u)
!     sym = guu(:,:,l)+transpose(guu(:,:,l)) ! symmetric gradient
!     do id=1,d
!       sym(id,id) = sym(id,id) + ldivu
!     enddo
!     fem(:,2:,l) = -rho(l)*nu(2,l)*sym
!
!     ! energy flux
!     fem(:,1 ,l) = -rho(l)*nu(1,l)*phc%cp*gtt(:,l) &
!                   + matmul( uu(:,l) , fem(:,2:,l) )
!   enddo
!
! contains
!
!  pure function tr(a)
!   real(wp), intent(in) :: a(:,:)
!   real(wp) :: tr
!   integer :: i
!
!    tr = a(1,1)
!    do i=2,size(a,1)
!      tr = tr + a(i,i)
!    enddo
!  end function tr
!
! end subroutine aniso_flux

!-----------------------------------------------------------------------

 pure function tr(a)
  real(wp), intent(in) :: a(:,:)
  real(wp) :: tr
  integer :: i

   tr = a(1,1)
   do i=2,minval(shape(a))
     tr = tr + a(i,i)
   enddo
 end function tr

!-----------------------------------------------------------------------

end module mod_turb_flux

